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We generalize the recently developed diagrammatic Monte Carlo techniques for quantum impurity 
models from an imaginary time to a Keldysh formalism suitable for real-time and nonequilibrium 
■ calculations. Both weak-coupling and strong-coupling based methods are introduced, analysed and 

' applied to the study of transport and relaxation dynamics in interacting quantum dots. 
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I. INTRODUCTION 



^NJ ' Quantum impurity models play a prominent role in nanoscience as mathematical representations of quantum dots, 
single-molecule devices and adatoms on surfaces. In general theoretical terms a quantum impurity model is a system 

with a finite-dimensional Hilbert space ("dot") coupled to one or more infinite systems ("baths") described by a 

I . Hilbert space with a continuum of energy levels. The equilibrium properties of quantum impurity models are by 
' now reasonably well understood theoretically and indeed in most cases the properties of interest can be computed 
, numerically to the necessary accuracy. 

' By contrast, the nonequilibrium properties of quantum impurity models are much less well understood. The subject 
QJ is of fundamental theoretical importance, as an instance of the basic problem of the properties of nonequilibrium 
■ quantum many-body systems. It is also of considerable experimental interest in connection with the properties of 
• ' quantum dots where the Kondo effect plays an important role in transport propertiesii^i^i^i^ Quantum impurity 
. models are also closely connected to the issue of transition rates and reaction dynamics in chemistry. 
' Quantum impurity models may be driven out of equilibrium in several ways. If a system is coupled to more than 
I . one reservoir, then a chemical potential or temperature difference between reservoirs can generate a nonequilibrium 
''O ' steady state in which current flows from one reservoir to another across the dot. One may also consider a transient 
^ ' or steady-state irradiation of the dot or the relaxation to steady state of an atypical initial condition; in either case 
one may have equilibrium or nonequilibrium reservoirs. While the basic formalism for dealing with these problems 
was established by Schwinger- and Keldysh^ in the early 1960s, and a wide variety of perturbative approaches have 
appeared (mainly tailored to specific physical applications), it is important to develop unbiased numerical methods 
which allow to test theoretical conjectures and to compare the properties of theoretical models to phenomena seen in 
experiments. 

Several numerical techniques have been applied to time dependent problems in interacting quantum dots. Numerical 
I renormalization group methods^ have been shown to provide impressively accurate treatments of relaxation dynamics 
fSj ' in dots with equilibrium baths and extensions to nonequilibrium baths have recently been proposed^^ However, 
— [ experience in equilibrium problems has been that these approaches, although powerful, are limited in the range of 
I . problems that can be treated and the range of energy scales that can be accessed. Path integral sampling techniques 
QQ ' introduced in the quantum chemistry context^ have recently been extended to the quantum dot problemiii These 
techniques require a finite "memory time" , and are therefore restricted to non-zero temperature and voltage bias. 
IL^ The time-dependent non-crossing approximation^^ gives access to long times and spectral functions, but is probably 
not reliable at strong interactions. The time-dependent density matrix renormalization group was also used to study 
the transport properties of quantum dots coupled to one-dimensional reservoirsji^iil 

In this paper, we present an extension to the nonequilibrium case of Monte Carlo approaches based on an unbiased 
sampling of diagrammatic expansion o^^d^d'^'d^ which, for equilibrium properties, have been shown to be powerful 
enough to access extremely low temperatures and flexible enough to treat a wide range of Hamiltonians. One of the 
specific implementations we present is closely related to recent work by Miihlbacher and Rabani^^ and Schiro and 
Fabrizic^*^ for non-interacting dots with coupling to phonons. An extension to interacting dots has been employed 
by Schmidt et al. in Rcf. [2l|. Here, we provide a systematic analysis of the real-time diagrammatic approach, 
including a discussion of the strengths and weaknesses of these methods, and the regimes in which accurate results 
can be obtained. We discuss the formalism and details of the measurement formulae and implementations, and present 
results for observables including dot double occupancy and current through the dot. The rest of this paper is organized 
as follows: in Section |TT] we outline the formalism we use and specify the model we treat, in Section IlIII we present 
the weak-coupling formalism and in Section IIVI the "strong coupling" or hybridization expansion method. Section IVl 
presents results for the time dependent dot occupation and double occupancy, and Section IVII discusses the current. 
Section IVIII is a conclusion and outlook. An Appendix presents derivations of some needed formulae. 
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FIG. 1: Example of a Monte Carlo configuration corresponding to perturbation order n = 5 and n+ = 3, n_ = 2. 

II. FORMALISM AND MODEL 
A. General considerations 

A quantum impurity model is described by the Hamiltonian 

Hqi — iJdot + ^^bath + -ffmix- (1) 

Here -ffdot describes a system with a finite dimensional Hilbert space, which we refer to as the "impurity" , or "dot" , 
^^bath describes one or more infinite reservoirs characterized by a continuum of levels, and i?mix the coupling between 
the impurity and the reservoirs. We assume that at time t = the state of the system is given by a density matrix po 
which will be specified in detail later. The statement that iJbath is an infinite reservoir implies that the distribution 
function describing the occupation of the energy levels of -ffbath is independent of the coupling to the dot. 
The theoretical task is to evaluate the expectation value {0{t)) of an operator O at time t, i.e. to compute 

(0(i)) ^Tr\poe'^o'^''"'^'^'"^Oe-'^o'it"HQ,{t")] ^2) 

(the generalization to operators with multiple time dependences is straightforward and will not be written explicitly). 
For a system in thermal equilibrium the issues in computing (O) are well understood. In this paper we are concerned 
with numerical approaches to describing the nonequilibrium situation. Nonequilibrium may enter through a time 
dependence of parameters in Hqj ("irradiation"), through the correlators of the operators in i?bath ("nonequilibrium 
reservoirs") or through an initial density matrix po which is different from the long-time limit. Our explicit considera- 
tions in this paper pertain mainly to the "nonequilibrium reservoirs" and "nonequilibrium po" case, but the methods 
we present generalize straightforwardly to other situations. 

One may^ view the expectation value in Eq. ^ as an evolution on the Schwinger-Keldysh contour illustrated in 
Fig. [T] from time t = (when the system is described by the density matrix po) to time t (at which the operator is 
measured), and then back to time 0. Our general strategy for evaluating Eq. ([2|) is to write Hqj as a sum of two terms: 
one, Hq for which the time evolution can be treated exactly and another, Hj, which is treated by a formal perturbative 
expansion. The expansion in Hj generates a series of diagrams which are sampled stochastically, using an importance 
sampling which accepts or rejects proposed diagrams on the basis of their contributions to (0) with, for example, 
= 1. Two forms of expansion are considered: One is a "weak coupling" method, in which H^ot is partitioned into 
a quadratic part H^^^^ and an interacting part Hjjj the combination iJ^^j -I- i?mix + -ffbath is diagonalized, po is taken 
to be the corresponding density matrix, and the expansion is constructed in terms of Hij. The other is a "strong 
coupling" (more properly, "hybridization") expansion in which -ffdot and -ffbath are treated exactly, po is the density 
matrix corresponding to the direct product of a lead density matrix and a density matrix describing the dot decoupled 
from the leads, and -ffmix is treated as a perturbation. The hybridization expansion for nonequilibrium problems was 
previously presented by Miihlbacher and Rabani^^ in the context of noninteracting electrons coupled to phonons, 
and has been applied to interacting dots in Ref. (21j . An essentially identical formalism has also been discussed in 
Ref. 0. 

Methods based on stochastically sampled diagrammatic expansions have had considerable success in equilibrium 
quantum impurity problems at temperature T > ,i5ii6,r7,i8 xhere, the expansion can be formulated on the imaginary 
time axis < t < 1/T (only one contour is needed) and the expansion parameter is —Hi{t) — e^^°{—Hi)e~'^^°. The 
fermionic sign problem can be avoided (at least for sufficiently small dots) and temperatures as low as 0.1% of the 
basic scales in the problem can be reached without inordinate effort. Three related sources of difficulty arise in the 
nonequilibrium problem. First, the expansion must be done for real times, so convergence of the perturbation theory 
is oscillatory rather than exponential: diagrams come with factors of i to powers relating to the perturbation order. 
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Second, two contours rather than one are required, doubhng the perturbation order required to reach a given time. 
Third, in nonequiUbrium situations the form of the density matrix is crucial to the quantities (such as the current) 
which are to be computed; thus it is essential that the computation proceed for long enough to build up the correct 
entanglement between the impurity and the bath. All of these factors limit the range over which accurate results can 
be obtained, but the crucial constraint is the dynamical sign problem resulting from the oscillatory convergence. 



B. Model 



The results in this paper are presented for the simplest possible situation, a "dot" consisting of a single spin- 
degenerate (a) level with a Hubbard interaction U, coupled by hybridization V to two reservoirs ("leads") labeled by 
a — L, R, with nonequilibrium entering via a possible difference between reservoir chemical potentials. The extension 
to more general situations is straightforward and involves no new conceptual issues. 

The Hamiltonian we consider is 

ifbath = ^ ^(e^,,-/ia)a^,t<., (3) 

a=L,R p,(T 

^-ix = Y.mK<^d„^h.c), (4) 

a=L,R p,iy 

Hit = {ed + U/2)J2nd,a, (5) 

(T 

Hu = U{nd^^nd,i- {nd,T +nd,i)/2). (6) 

It is also convenient to define 



Hdot - -f^dot + Hu- (7) 

The initial density matrix is such that the correlators of lead operators are {frix) = {e^^^ + 1)^^ is the Fermi 
distribution function for temperature T) 

and the statement that ffbath describes infinite reservoirs is the statement that Eq. ^ holds at all times. 

The model has three important energy scales: which controls the steady state dot occupancy, the interaction 
scale U, and the level broadening 

r"(c.) = 7r5]|^;^p5(c.-e^) (9) 

p 

associated with lead a. The total level broadening is 

r + (lo) 

and the dimensionless measure of interaction strength is U /T. Very roughly, strong coupling physics appears for 
[/ > ttF while the opposite limit is reasonably well described by perturbation theory in U (see Section lv|). 



III. WEAK-COUPLING ALGORITHM 



A. Weak-coupling expansion and auxiliary field decomposition 



In the weak coupling expansion we treat Hq = + i?mix + ^^bath exactly and Hjj as a perturbation. Hq is a 
noninteracting problem for which the density matrix and all correlators of the dot-lead system can be determined 
exactly. We take the initial density matrix to be the steady-state density matrix corresponding to Hq (here we assume 
the temperatures of the two leads are identical; the generalization to unequal temperatures is straightforward) 
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and consider the interaction to be turned on at time t = 0. 

We formulate the perturbation theory in J7 as a real-time incarnation of the recently developed continuous-time 
auxiliary field method of Ref . [l^ , which itself is an adaptation of ideas in Ref. [1^ to impurity models. The starting 
point for the real-time auxiliary field method is the following expression for the identity: 

1 = Trpoe''^^"°+""^'^/^^e''^^"°^""^^/^\ (12) 

with K a constant which is in principle arbitrary and may be chosen to optimize the simulation. As discussed below we 
find that choosing K to be negative, and small in magnitude appears to work best. Using an interaction representation 
in which the time evolution of the operators is given by 0{s) = qIsHoq^-isHq ^^^^ rewrite Eq. (fT^ as 



1 = 



with T the time ordering and T the anti-time ordering operator, and expand the time ordered exponentials into a 
power series. This leads to the expression 

1 = Trp^Y^{-iK/tY' [ dii...[ dt„e**i^"(l-iiJa/-?f)---e'(*"~*"-^^^"(l-i-ffc//^)e'^*"*"^^° 

X Y^{iK/t)'^ f dti... f dt„e-*(*-*")^° (1 - tHu/K) . . . e-*(*^-*i)^° (1 - tHu / K)e-'*'"° . (14) 

Using the explicit form for Hu (Eq. ([6])) and the auxiliary field decomposition of Ref. [1^ we can rewrite the interaction 
term as 

l-{tU/K)in,^^n,,^-{n,,^+n,,^)/2)) - 1/2 ^ e^^("'''T-".a), (15) 

cosh(7) = l + {tU)/{2K). (16) 

Note that the constant K has been introduced to enable this decomposition. The trace is now a product of exponentials 
of one-body operators, 

1 = ^^(-*)"V(if/2<)'"+" ^ fdh-.-t di„j'dh...f dU^Hil/Tre-^''"-) 

X Tr gPI^'^'^"'^-'' gi(tm-tm-l)-ffo,<Tg7SmCrrad,„g-i(tm-t„)/fo,<Tg7Sn<^"d,<T g" « (*2 " * 1 ) -^0 , (t 1 (^"d , <t g 1 -H" , ,t 

and can be expressed^^ in terms of determinants of two (n -I- m) x (n -\- m) matrices 

^ e^' - iGoA^^' ~ I) (17) 



as 



1 = ^^(-z)'"z"(if/2i)"+" ^ J2 dtmfdh...f dt,,YldetN-\ (18) 

™ n 5i,...,5„si,...,s„.-^0 Jt^-, Jo Jt„-, 

with e^- = diag(e'^^i'", . . . , e^^"'", e*"^" , e^^^'"), and Go,a given by 

r d' i" \ — i ^o,(t(^'' ^")' ^'k < ^'k /iq\ 
^oAiK.tK) - I G>^(t',i"), t'x > t'ii ■ ^^^^ 

In the above expression, G^{t,t') = i{d'<{t')d{t))o, G^{t,t') = -i{d{t)d^{t'))o, Ik is the "Keldysh time" coordinate 
along the unfolded Keldysh contour (Fig. [T|) and t the time corresponding to Ik- These Green's functions may be 
computed by standard methods. A general expression is presented in Appendix A; for our actual computations we 
will use the infinite bandwidth limit in which the level broadening is independent of w so that 

G</>rt' t") - ±i V F" / '^^c-»^f*'-«"1 iTtanh(^) 

Go {t,t) - ±1 I / (c. - - [7/2)2 + F2 ^^^> 



a=L,R 



with the upper sign pertaining to Gq and the lower sign to G^ 
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B. Detailed balance and fast updates 

The algorithm samples auxiliary Ising spin configurations {(tx.!, si), {tK,2, S2), ■ ■ ■ {tK.n, Sn)} time ordered along 
the "Keldysh" contour ^ i — > (see Fig.[T]) by random insertions and removals of spins. The complex "weight" of 
a spin configuration is given by 

w{{{tK,l, Si), {tK^2, S2), . . . {tK,n,Sn)}) = (-i^^ ) (i"+ ) (i^di/2t)"- [] ' (21) 

a 

where n+ denotes the number of spins on the forward contour and n_ the number of spins on the backward contour 
{n = + n_). 

The detailed balance condition for insertion/removal of a spin is similar to the imaginary time formulation of 
Ref. [l8| . Assuming that we pick a random time on the unfolded contour of length 2t and a random direction for this 
new spin (pP''°P(ri — 1 ^ n) = {\/2){dt/ (2/;)), and propose to remove this spin with probability pP''°P(n n — 1) = l/n 
we get 

p{n~l^n) _ ^2Ky. det(jV-i)^ 

P(n-n-l) ~ n lidet(7V-_\).' ^''^ 

with the factor +z corresponding to a spin which is inserted on the forward contour and — i to a spin which is inserted 
on the backward contour. 

For the fast updates, let us consider the most complicated case, which is the insertion of a spin. This update adds 
one row and one column to the (n — 1) x (n — 1) matrix N , resulting in the n x n matrix N' (we assume here that 
this new row/column is the last one, n, and drop the spin index). The determinant ratio is 



det(7V'-i) 



det(7V-i) 



^ = (e^ - iGoie^ - I))n,n - J2 ^^(^^ - - ^)kn, (23) 



i=l 



with Ri = (^^ ~ iGo{e^ — I))njNj^i. The calculation of this quantity requires 0{n^) operations. The new 
matrix elements are given by 

Nl^=N,^j + ^L,R,, (24) 

Kn = -^L,, (25) 

7V;, (26) 

Kn - (27) 

r 

with i = 1, . . . , n - 1 and Li = J^'jZi ^ij(e'^ - iGo{e^ - I))j,n- 

From Eq. ([27]) it follows that computing the determinant ratio for removing a spin is 0(1). The elements of the 
reduced matrix are obtained as 

N'- N' ■ 

N.j-K-^§r^- (28) 

C. Green's function, dot population and double occupancy 

To measure the Green's function Ga{t'j(,tj^) we have to insert an operator at time t'^^ and an operator dj. at 
time t'l^. The weights of these configurations w{{{tK.i, si), . . . {tK,n,s„)}; d^{t'j^)d^„{t'l^)) are related to those defined 
in Eq. ([21]) by 

(29) 



w{{{tK.l,Sl),...{tK,n,Sn)};d{t')d\t")) _ 1 



j{{{tK,l, Si), . . . {tK,n, Sn)}) det Na 



-1 





iGo,cr{tK,i,tx) ' 


-iGoAtK^tK^jKe-^"'^ 


-1) 


iGo,cr{t'K,tK) , 
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Hence, the Green's function can be obtained as the Monte Carlo average of the quantity (see also Ref. [18']) 



i,3 = l 

which yields the measurement formulas 

n^n^itK) = {{l-iG^{tK,tK)Kl~iG^{tK,tK))). 



(30) 



(31) 
(32) 
(33) 



D. Current measurement 

The current from the dot to the left lead is 



(34) 



Thus, in terms of the composite lead operator aj^ ^ = X^pei ^pa^^\i ^'^'^ 

= -2Im V V(-*)"z"(i^/2t)™+" V V /* dh ... t di„, f dh... f dtn 



(35) 



with a the spin which is opposite to a (this spin component has no operator a and thus simply gives the usual factor 
detN^^). The measurement of the current is thus very similar to the measurement of the Green's functions, but one 
factor in the Wick decomposition is now 



^ 1 A>(t,t') = -{dAt)&^^{t'))„, tK > I'k ' 



(36) 



A derivation and a general expression are given in Appendix A. In the infinite bandwidth limit we have 

dio 



A>{t,t') 



2tt 



2tt 



f{uj - hl) 



{uj-ed-U/2r + r^ I (/(c.-ml)-1) 



(37) 



The trace factor in Eq. (I35p for an 7i-th order diagram corresponding to the n x n matrix ^ is the determinant 
of the {n + 1) X (?i + 1) matrix 





A{tK,r,t) 


-iGoAt,tK,j){e-"''^ 


-1) 


Ait,t) 



The current can thus be expressed as follows: 



II = -2Im^^u;f- = -2Im^ 



T,c\'Wc\iwi' /\Wc\) 



Wc\ I \wa\ ('/>c)|u;,| 



(38) 



(39) 
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with (t)c the phase of the weight Wc (Eq. (|2T|) ') and 



wi' det/^^MetM, i ^ A{t,t) +y^iG^At.tK,nW' - l)NX^„A{tK,n.,t). (40) 
detiV^ detiVCT 



Combining Eqs. (j39p and (|40p . the current measurement formula becomes 

/ = =-2Imy [A(t,i) + (ViGo,.(t,tif,„)[(e^- -l)iV,]„,™A(iK,m,0'^c) 7-^ 1- (41) 

The first term in this expression is the steady-state current for the non-interacting system 

/o = -2Im(2A(M)) = Sy (^_,^_t//2)2 + r^ ■ (42) 

E. Real-time Hirsch-Fye method 

In order to assess the efficiency of the continuous-time weak-couphng approach, we have also performed calculations 
using the real-time version of the Hirsch-Fye method.^** In this method, time is discretized along the Schwinger-Keldysh 
contour and the identity is expressed as 

L/2 L/2 

'[Ho+Hu] ^43^ 

(44) 

where we used the Trotter breakup and L denotes the (even) number of time slices. In the real time Hirsch-Fye 
method, the interaction term for each time slice is decoupled using the Hubbard-Stratonovich transformation 





L/2 


L/2 


= Trpoe'*^e-**^ = 


: Trpo n e 


At[Ho+Hu] -Q g 
(=1 


L/2 


L/2 




1=1 


AtHu -Q g- 
1=1 


-iAtHa -iAtHu 



" ^ - e^''("T-"i)^ A = cosh-^e*'^*^/^). (45) 



After this decoupling, we obtain a product of exponentials of one-body operators and the trace can thus be computed 
analytically. Besides the fact that A is a complex number and that the non-interacting Green's function is given by 
Eq. (|19p . the derivation of the algorithm and the sampling procedure are identical to the original imaginary-time 
Hirsch-Fye method. Two sources of error exist in this method. One is the discretization error due to the Trotter 
breakup and the other is the stochastic error which becomes severe at large L due to the sign problem. Because of 
these limitations, the real-time Hirsch-Fye method is restricted to short time calculations, and indeed the time limits 
appear to be more stringent than in the continuous-time methods we introduce here (see for example Fig. [8]). 

IV. HYBRIDIZATION-EXPANSION ALGORITHM 
A. Formalism 

A complementary diagrammatic Monte Carlo algorithm can be obtained by performing an expansion in powers of 
the dot-lead hybridizations V. This simulation approach has been introduced for equilibrium systems (imaginary- 
time formalism) in Refs. [T6llT7ll2^ and was recently discussed for a noncquilibrium dot with phonons (but without 
electron-electron interactions) in Refs. [l9ll20l |. It has been applied to interacting dots in Ref. |2l|. We will present 
here the derivation for the impurity model defined in Eqs. (O-®, but the method can easily be extended to general 
classes of impurity models by using the matrix formulation of Ref. 17]. 

In the hybridization expansion approach one adopts an interaction representation with respect to the dot-lead 
mixing, so the time evolution of the operators is given by the local part of the Hamiltonian, iJioc = -ffdot + ^^bath, and 
the starting point is the identity 



1 = 
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The initial state of the system is specified by the density matrix po = Pdot 'Si Pbath, with pbath a function of inverse 
temperature P and the chemical potentials Hl.r- the calculations presented here we assume that the dot is initially 
empty, pi,„p = |0)(0|. 

Expanding the time ordered exponentials into a power series yields 

1 = Trpo J2 / dii-- - [ diraH^Uii) • ■ ■ i?mix(im) 
Jo Jt„,-i 



(47) 



Because = E.iH^i. + <l) with iJ,^^, = Ec.=L..REpV,^a^},d,, <i = and the time evolution 

conserves the spin, we need for each a separately an equal number of creation and annihilation operators on the 
Keldysh contour ^ t ^ 0: 



1= E 

m„+ri„=m' +Ti' (J 



X 



dfL 



(-0 



dt'^ 



dC. / dti 



dtl 



dt'^ 



dt 



X Tr 



(48) 



where T is the anti-time ordering operator for the ts and T the time ordering operator for the ts. At this stage we 

rs a'^ ,^ i 



can separate the bath operators ^ from the dot operators da and write 



1 



m„+n„=m' +n' a 



X 



dt'-, / dt\ 



dtl^ dt'l 



dt'Z 



Pdot 



f rn d.(t?)4(iT)rf.(i2 )4(ir ) • • ■ e^^-*e-^-* . . . d.(t^)4(t-)d.(t?)4(tr) 



X Tri 



bath 



p.-fTn E E E E 



Pi p\ • ■ ■ ^Pl "^p'j 



at(tl)a,(f-)aT(t-)a.(tr)...e 



(49) 



with ai G Since the leads are non-interacting (see Eq. ^) we can evaluate the factor Trbath[- ■ •] exactly. 

Due to Wick's theorem one obtains a product of two determinants det M~^, with the size of given by the 
number of operators d^ on the Keldysh contour (toct + n^). The matrix elements are given 

M-\t,j)=tA{ej,^,,t'^^^), (50) 

where t J- ^ denotes the position of the ith annihilation operator and t'^ j the position of the jth creation operator for 
spin a on the unfolded Keldysh contour. The function A is discussed in Appendix A and is 



with 



^ \ A<{f -t)^ A<{t' -t) + A<{t' ~t) tK> t'^, 



/CO 1 
-^^e~-*r«(.;)/(^-^„ 

A>(t) 



2i 



*^"(c.)(l-/(a;-Ma))• 



(51) 

(52) 
(53) 
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FIG. 2: Illustration of the band cutoffs considered in the simulations. The soft cutoff (left panel) is an exponentially decaying 
r{uj) which is centered at the chemical potential. The band with hard cutoff (right panel) is symmetric around lo = and does 
not shift with the chemical potential. 



We give here simple expressions for two types of bands, illustrated in Fig. ^ which are exact in the limit T — > 
and a good approximation for T iUc- The first is a band with soft cutoff, centered on the chemical potential, 



rs,ft(u;) = r"e-i-^°i/-=. 

For symmetric voltage bias {fiL — ^fJ-R — ^/2) and symmetric couplings (F^ — F/j) we obtain 



cos(-jt) 



The second example is a flat band centered at zero with a hard (Fermi- function like) cutoff a,t uj — ±Wc, 

rhard(^) 

which yields 



pa 



(1 



-))(! + e-'^("+"<=))' 



^haid 



it) 



F 



cos(l^t) 
/3sinh(|i) ~ j/sinh(fi) 



(54) 



(55) 



(56) 



(57) 



To evaluate the trace over the impurity states in Eq. (|49p . Tr^l- . . ], it is useful to employ the segment representation 
introduced for impurity models with density density interactions in Ref. p^ . The sequence of dot creation and 
annihilation operators uniquely determines the occupation of the dot at each time, and we can represent the time 
evolution using collections of segments for spin up and down electrons as shown in Fig. [31 Each segment depicts a 
time interval for which an electron with corresponding spin resides on the dot. The trace over the impurity states can 
then simply be expressed as 



= Pimp(c) exp 



/ J V^forward ^backward 



overlap 
forward 



T overlap 
backward 



(58) 



Here, Pimp(c) is the element of the impurity density matrix which is compatible with the operator sequence c = 
\t%^ . . . , m^+n„ It ■ ■ ■ m' +n' } (assumed here to be 1 for configurations which start and end with an empty 
dot and zero otherwise), l'^ the length of the segments for spin cr and ;°™riap length of the overlap between spin 
up and down segments. 

Hence, the Monte Carlo simulation samples collections c of segments on the doubled "Keldysh contour" (one for 
each spin) according to their weight 



w{c) = JJi™'+™-(-i)"-+"-detM-idf 



X Pimp (c) exp 



forward 



backward/ 



Ty / 70verlap lovcrlap 
V forward backward'' 



(59) 
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FIG. 3: Segment representation of a Monte Carlo configuration corresponding to perturbation order 4 for spin up (upper 
contour) and 2 for spin down (lower contour). Dot creation operators are shown as full circles and annihilation operators as 
open circles. The segments represent the time intervals in which an electron of the corresponding spin resides on the dot. 




FIG. 4: Segment configurations obtained from the expansion of the current (Eq. (|62|) ') in powers of the dot-lead hybridization. 
There is a fixed operator do- (red open circle) at time t and the hybridization functions connecting to this operator have only 
a left (i) component. 



We implemented the following local updates of the segment configurations: i) insertion/removal of a segment, ii) 
insertion/removal of an anti-segment (empty space between segments) and iii) shifts of segment end-points. In order 
to use fast update formulas similar to those discussed in Section IlIII we store and manipulate the matrices Ma, that 
is, the inverse of the matrices defined in Eq. ([50]) . 



B. Measurement of the Green's function, density and double occupancy 



The Green's functions can be obtained from the matrix M in a procedure analogous to the one proposed for 
imaginary-time simulations in Ref. [l6j . Particularly simple is the calculation of the density and double occupancy. 
From the segment representation it immediately follows that Uait) is the probability to have a segment of spin a 
present at time t, while n|ri|(i) is the probability to find overlapping segments at time t (taking into account the 
signs of the Monte Carlo configurations): 

((/«c(5(segment of type cr at i))|,„^l 
na{t) = j-T (60) 

((/)c'5(segments of type T and J, at t))|^^| 
n^niit) = — (61) 



C. Current measurement 



The current Jlct = ^Slm^^g^ Vj^^lapJ^dc) = — 2Im(a|^ ^.d) can be measured as explained in Ref. [T3|. We expand 
the quantity 

lL,{t) = -2ImTrpo(Te*^o (^))e^*^i"a[^d,e-^*^i"(re-^^*'^^^— (62) 

in powers of -ffmix, which leads to the same collection of diagrams as discussed above, except that there is now 
an operator d^r fixed at time t and that the hybridization functions A connecting to this operator have only an 
i-component: 

M-i(i,j) = zAi(if^^„i^,^.)+iAfi(i^^„i^,,)(l-^t,t^,J- (63) 
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Having identified the Monte Carlo configurations c (illustrated in Fig. |4]) and their weights Wc we can implement a 
random walk based on \wc\ and measure the current as 

iLa =^Wc^ {<l>c)\w,\^\Wc\- (64) 
c c 

In contrast to the density measurement (which was based on an expansion of the identity so that = 1/ (0c)|uic|) 

we cannot directly measure the normalization factor \wc\- One possibility to get rid of this unknown factor is 

to consider the ratio I/I^^'> between the current and the lowest order contribution I^^^ which can be calculated 
analytically. Since 

IlI = {<l>cS{c 1st order))|^j ^ \w,\ (65) 

c 

we can measure the current as 

^-(0,(5(clstorder))|„^|- ^""^ 



V. RESULTS: PERTURBATION ORDER, DENSITY AND DOUBLE OCCUPANCY 



A. Perturbation order and average sign 



The average perturbation order increases linearly with the time interval to be simulated, and is not per se the 
important limiting factor in the simulations. The main constraint is a dynamical sign problem: the factors of 
(±i) associated with each order of the expansion and the complex determinants mean that the average sign of the 
diagrams contributing to any quantity decays exponentially as the perturbation order is increased . These phenomena 
are illustrated in Fig. [5] which presents results obtained using the hybridization expansion algorithm on a model of 
spinless fermions. The same behavior is found in interacting models and in the weak coupling algorithm. The left panel 
shows the distribution of perturbation orders for simulations over different time intervals. The mean perturbation 
order can be estimated from the positions of the maxima in these curves. The right panel shows the average sign, 
which is seen to decay exponentially with the length of the time interval to be simulated. Note that both diagrammatic 
algorithms can treat temperature T — without particular difficulties. 

Accurate measurements of physical quantities can be obtained for (sign) > 0.001, and whether steady state can be 
reached depends on the method, the parameters, and the observable. Non-zero temperature and voltage bias tend 
to reduce the sign problem, but not enough to enable simulations on significantly longer contours. The important 
effect of a non-vanishing voltage bias is to accelerate the convergence to steady state, at least in the weak-coupling 
approach. As can be seen e.g. from the right hand panel of Fig.[5l while the basic scaling behavior is the same for all 
models and parameters, prefactors can depend substantially on details. A careful effort to optimize parameters has 
not yet been undertaken, but seems likely to be worthwhile. 

The left panel of Fig. [S] shows that in the weak-coupling approach, the average perturbation order (at fixed t) 
depends on the interaction strength. As in the imaginary-time version of this algorithm,— the perturbation order 
grows roughly linearly with increasing [/, making it difficult to study dots with U/T > 3. Larger values of K also 
lead to a larger perturbation order and hence to a more severe sign problem, as illustrated in Fig. [T] The parameter 
K can also be chosen negative or complex (we used K = —0.01 in all our weak-coupling simulations). In this case, 
7 is complex and some of the phase oscillations are shifted from the (±zi(r)"± to the determinant. It is in principle 
possible to choose different constants K+ and if_ on the forward and backward contour. If = iK = —K^ (with 
K positive), then all the phase oscillations from the {±iK±)'^^ are eliminated. However, it turns out that — 
leads to a perturbation order which is about the same as for iK on both the forward and backward contour, which 
in turn is somewhat worse than ~K on both the forward and backward contour, so that there appears to be no 
particular advantage in this choice of K. 

The right panel of Fig. [S] shows that in the hybridization expansion algorithm, the average perturbation order is 
essentially independent of interaction strength. This is in contrast to the imaginary-time version of this algorithm^ 
where the perturbation order decreases with increasing interaction strength. From Eq. (j59p it follows that the 
interaction term merely adds a phase to the Monte Carlo weight and therefore does not affect |w(c)|. While the 
algorithm can treat strong interactions, it is limited to finite bandwidth, since the average perturbation order diverges 
as the bandwidth goes to infinity (this is the reason for the dependence on cutoff seen in the right hand panel of 
Fig. [5]). We find, however, that systems with larger cutoff reach steady state more rapidly (as in Fig. 1 of Ref. 21]). 
We have not yet attempted to optimize the cutoff to strike the best balance between perturbation order and time 
needed to reach steady state. Such an optimization would be worth while. 
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FIG. 5: Distribution of perturbation orders and average sign obtained using the hybridization expansion algorithm for a non- 
interacting dot with soft cutoff, T — 0, V = 0, ea/T = —0.5 and a single species of fermion. Left panel: distribution of 
perturbation orders for different lengths of the contour {tV = 1.25, 1.50, ... ,3 from left to right) and cutoff Wc/F = 40. The 
average perturbation order grows ~ t. Right panel: average sign as a function of time for indicated values of the cutoff. 




perturbation order perturbation order 

FIG. 6: Distribution of perturbation orders for different values of the interaction U,T — 0,V — 0, ed + U/2 — 0. The left panel 
shows the distribution of perturbation orders for the weak-coupling algorithm {tT — 1, infinite bandwidth), where the average 
perturbation order grows ~ U. Right panel: distribution of perturbation orders for the hybridization expansion algorithm 
{tV = 1.5, ojc/F = 10). Here there is almost no dependence on interaction strength. 



B. Density and double occupancy 

The left panel of Fig. [8] shows as black symbols the evolution of the double occupancy obtained using the weak- 
coupling algorithm in equilibrium {V = 0) and at zero temperature for a system tuned to be at half filling. At time 
t = the double occupancy takes the value 0.25 appropriate to the noninteracting half filled system. The effect of 
the interactions is to reduce it. The right panel shows the dot occupancy per spin, computed for a level position 
corresponding to a quarter- filled dot at U ~ 0. Turning on the interaction increases the dot occupancy; this is a 
precursor of the Coulomb blockade plateau. 

One sees from the figure that for U/T = 2 it is possible to obtain a good estimate of the steady state value, whereas 
for U/r = A the perturbation order grows too rapidly with t and the sign problem becomes severe before the system 
approaches the steady state. Whether or not steady state can be reached depends on the observable. The statistics 
for the density is significantly better than for the double occupancy, so density calculations can be carried to longer 
times. 

We also show in Fig. [5] results obtained with the real-time Hirsch-Fye method for tT = 1 and different numbers 
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FIG. 7: Left panel: perturbation order distribution for tV = 1, U/T — 2, V/T = and several positive values of K. The average 
signs in these simulations are 0.24 {K = 0.1), 0.04 {K = 1) and 0.0002 (K = 10). Right panel: complex K of norm 0.01. The 
best choice appears to be the negative K. 
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FIG. 8: (color online) Weak-coupling results for the double occupancy and density computed at T = and = with 
parameter K = —0.01 and for infinite bandwidth. Left panel: double occupancy for U/T — 2, 3, 4 at half filling (e^ + J7/2 = 0). 
Right panel: density per spin for ed + U /2 — T with the other parameters the same. The full dots show the results from the 
real-time Hirsch-Fye method for tV = 1 and indicated values of the number of time slices L. 



of time slices. As the number of time slices is increased, the systematic error due to the Trotter break-up decreases 
and the result approaches the continuous-time curves (which are free of systematic errors) . Comparison of the left 
and right panel shows that the density is less sensitive to Trotter errors than the double-occupancy. Because the 
sign problem in the real-time Hirsch-Fye method becomes severe for L > 30, longer times can only be reached at the 
expense of larger discretization errors. In our calculations we found that the weak-coupling continuous-time algorithm 
allows to roughly double the time interval which can be simulated, compared to Hirsch-Fye. 

Simulations a,t V = and T — suffer from the most severe sign problem. At non-zero voltage bias, the system 
reaches steady state more rapidly, as illustrated in the left hand panel of Fig. [51 For U/T — 2 and 3 we can therefore 
obtain an accurate estimate of the steady state double occupancy. The voltage dependence of this quantity is plotted 
in the right hand panel of Fig. [HI As the voltage bias is increased, the steady state double occupancy drops, reaches 
a minimum and then increases with increasing V toward the non-interacting value of 0.25. The initial drop in the 
double occupancy is the result of the destruction of Fermi liquid coherence with increasing voltage bias, similar to 
the destruction caused in equilibrium by a non-zero temperature. At larger voltage bias a reversion towards the non- 
interacting value of 0.25 is evident. This non-monotonic behavior was also observed by the time-dependent density 
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FIG. 9: Voltage dependence of the double occupancy for an infinite flat band, T = 0, half filling {ed + U/2 = 0). The left panel 
shows the time evolution of the double occupancy for U/T = 2 and indicated values of V. At finite voltage bias, the system 
reaches steady state much more rapidly than for V = 0. Right panel: steady state value of the double occupancy as a function 
of voltage bias for U/F — 2 and 3, T = and half filling. 
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FIG. 10: Hybridization expansion results for the double occupancy and density. Left panel: double occupancy for indicated 
values of U, band cutoff ojc/F — 10, f3T = vT = 10, V = G, ed + U /2 = 0. Right panel: density per spin for the same parameters. 
The initial state is an empty dot decoupled from the leads and at f = the dot-lead hybridization is turned on. 



matrix renormalization group methodic 

Figure [TO] shows results for double occupancy and dot occupation obtained using the hybridization expansion 
algorithm. Here, the initial state is an empty dot which is decoupled from the leads and at i = we turn on the 
hybridization. The left panel shows the time evolution of the double occupancy in a dot with + U/2 = 0, a hard 
cutoff LOc/T — 10, /3r — vT ~ 10, and the right panel shows the evolution of the density per spin. Increasing the 
interaction accelerates the approach to equilibrium, and leads to a slight overshooting of n[t). The reduction of the 
double occupancy is roughly consistent with the result from the weak-coupling simulation (note that the band widths 
are different). 
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FIG. 11: Current computed as a function of voltage bias in the infinite bandwidth model using 4*'' order perturbation theory 
in the interaction U (from Ref. [30l]l. 



VI. RESULTS: CURRENT 
A. Qualitative Picture: perturbative and mean field results 

To orient the discussion of our results for the current we present here a brief outUnc of the expected qualitative 
behavior, along with perturbative and mean-field calculations. In the noninteracting limit, the d-density of states of 
the model defined by Eq. ([T]) takes an approximately Lorentzian form with a peak at the d-level energy Sd and a 
width of order F. As U is increased, the structure of the d-density of states changes: the peak broadens, and at large 
enough U splits into two. The density of states in the region between the two peaks becomes small, except that in 
the ^ 0, T — > limit a narrow peak (the Kondo resonance) appears at the Fermi level, so the Fermi level density 
of states remains essentially unrenormalized. At temperatures or voltage biases greater than the Kondo scale (which 
becomes exponentially small at strong couplings) the Kondo peak is believed to be destroyed, leaving only the small 
density of states (Coulomb blockade) behavior. 

The current / is, up to various constants, given by the integral over the voltage window ~V/2 < e < V/2 of the 
product of FiFjj/F^ and the density of states. In the noninteracting case / starts out linearly with V and saturates 
for V ^ r. As [/ is increased the broadening of the peak means that the y-value needed to reach current saturation 
increases. The Kondo physics implies that the T = linear response current is essentially independent of interactions, 
but what happens at larger V in the strongly correlated regime is unclear. 

Figure [TT] shows the current computed from 4*'* order perturbation theory in U by Fujii and Ueda'^'' for the infinite 
bandwidth version of Eq. ([T]). The initial linear rise and eventual saturation of the current are clearly visible, as is 
the increase in the saturation voltage as U is increased. Also visible in the calculation is the i7-independence of the 
linear-response current and hints of the formation of the Coulomb blockade plateau at intermediate V and larger U. 
Of course, the reliability of low-order perturbation theory at these interaction strengths may be questioned. 

Figure [12] shows the results of computations performed using mean field theory^^ as well as phenomenological 
generalizations. Details of the calculations are given in the Appendix, but the essence is as follows. In mean field 
theory of the model studied here, at T = and in equilibrium, a transition occurs at Uc = ttF between an unpolarized 
weak coupling state and a strong coupling state characterized by a frozen local moment and a spectral function split 
into upper and lower Hubbard bands. This is the mean field representation of Coulomb blockade. The mean field 
theory does not capture the Kondo effect, so for U > Uc the near Fermi surface density of states is simply suppressed. 
For U > Uc, as the voltage is increased, the degree of spin polarization decreases and within mean field theory we find 
a sharp phase transition at which the properties revert to those of the unpolarized state. The dashed-dotted curve 
in the inset of Fig. [T^] shows that for the model studied the transition is first order (jump in to) and for U = 12F 
occurs a.t V ~ 7F. Ref. [11] presented analytical arguments that a polarized phase would extend to infinite voltage, 
at least in a model with an infinite bandwidth; this behavior is not found in our numerical solution of the finite- 
bandwidth model. The main panel of Fig. [12] shows the mean field current computed in various approximations. 
The dashed double-dotted trace (black on-line) shows the current a.t U — (the small differences from the U = 
trace in Fig. [TT] arise because in Fig. [T^ a finite bandwidth is used, whereas in Fig. [TT] an infinite-bandwidth limit 
is taken). The dashed-dotted curve (red on-line) shows the current computed from mean field theory. Comparison 
to the noninteracting (dashed double-dotted) curve reveals the Coulomb-blockade suppression of the current at small 
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FIG. 12: Main panel: dashed-dotted line (red on-line): current computed in mean field theory for U — 12T as a function 
of voltage bias V assuming negligible pseudothermal broadening. Dotted line (red on-line): current computed in mean field 
theory assuming a pseudothermal broadening equal to 20% of the voltage bias. Solid line (blue on-line): current computed 
using mean field theory with gap fixed at the V = value and negligible pseudothermal broadening; dashed line (blue on- 
line): current computed using mean field theory with gap fixed at the V — value and pseudothermal broadening equal to 
20% of the voltage bias. Dash double-dotted line (black on-line): current computed for the non-interacting model (negligible 
pseudothermal broadening). Inset: dot magnetization as a function of voltage bias V/F computed in mean field theory for 
negligible pseudothermal broadening (dot-dashed line, red-on line) and moderate pseudothermal broadening (dotted line, red 
on-line). All computations were performed for a hard cutoff with ljc = lOF and i^F = 10. 



bias and the reversion to the noninteracting result at higher bias. Within mean field theory the reversion occurs via 
a first order transition at a critical bias of the order of one half of the Coulomb gap. Mean field theory is of course 
not an entirely accurate description. For example, as noted by the authors of Ref. [28], the transition is an artifact of 
mean field theory. One would expect features of the Coulomb gap to persist at high voltage biases. To qualitatively 
assess the consequences of this physics we show as the solid line (blue on-line) the results of a computation in which 
the Coulomb gap (splitting of the density of states into two peaks) is fixed at its F = value. A broader range of 
current suppression and a high saturation voltage are evident. 

The mean field theory is deficient in an additional way. Ref. [1^ showed that mean field theory misses the fact 
that bias voltage functions as an effective temperature (proportional to the bias voltage times a numerical factor 
related to scattering phase shifts) which broadens all of the properties. In order to asses the qualitative effect of this 
consideration we modeled the pseudothermal broadening effect of a voltage bias by performing the calculations at a 
temperature chosen to be T^g = 0.2V. The dotted lines in the inset and main panel show that the pseudothermal 
broadening effect converts the first order transition into a second order one. More significantly, we see that including 
an effective temperature tends to decrease the current at higher biases. The numerical calculations discussed below 
will be seen to be most consistent with the fixed gap, pseudothermally broadened mean field calculations. 



B. QMC results: Current 

1. Weak Coupling Expansion 

In the weak-coupling simulations, the current starts from the steady-state value for the non-interacting dot, 
— 41mA(0,0), decreases in magnitude after the interactions are turned on, and eventually converges at sufficiently 
long times at the value corresponding to the steady state current through the interacting dot. As shown in Fig. 1131 
useful estimates for this steady state value can be obtained for interaction strengths U/T < 3. While the first term 
in Eq. ([57)) can be computed directly for the wide band limit, the second integral requires a frequency cutoff ujc- 
However, we found that the current results are insensitive to this cutoff-value, as long as Uc ^ V. All our results were 
obtained for ujc/T = 10. 

Figure [T3] presents the time dependence of the current obtained from the weak coupling algorithm for different 
values of the voltage bias at U = 2T (left panel) and at different interaction strengths for fixed bias V = 4T (right 
panel). We see that the effect of the interaction is to reduce the magnitude of the current. However, the corrections 
are relatively small at the U and bias voltages studied: the noninteracting systems already gives a good approximation 
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FIG. 13: Weak-coupling expansion results for the current at temperature T = 0. The left panel shows the time evolution 
for U/r = 2 and indicated values of the voltage bias. The value at time t — corresponds to the non-interacting current 
lo = — 4Imyl(0,0). As the interaction is turned on, the magnitude of the current decreases and eventually converges to the 
steady state value for the interacting dot. The right panel shows the time dependence of the current for V/T = 4 and U/T = 2 
and 3. 
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FIG. 14: Weak coupling expansion results for the voltage dependence of the steady state current for T — and indicated 
values of the interaction U. The left panel shows the total current and the right panel the reduction in the steady state 
current produced by the interaction, I{U) — /(O). The interaction correction peaks at a value of V which is comparable to 
the interaction, or somewhat larger than the twice the "level broadening" F. The thick black line in the right panel shows the 
prediction for U/T — 2 from 4**^ order perturbation theory.— 



to the current if the interaction is not too strong. 

Figure [M] plots the current as a function of voltage bias for different values of U. As shown in the right hand 
panel, the largest reduction of the current is observed for V/T « 2.5, which is comparable to the interaction strengths 
U/T = 2, 3. As a consistency check, we show in the right hand panel as a thick black curve the interaction correction 
for U/T = 2 deduced from the 4*'' order perturbation calculation of Ref. [30|. The perfect agreement with the Monte 
Carlo data shows that the perturbative calculation gives accurate results at this small coupling strength. 



2. Hybridization Expansion 



In the hybridization expansion method (as it has been implemented here) the initial state is an empty dot decoupled 
from the leads so that at i = there is no current. As time evolves from t = the current must build up to its steady 
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FIG. 15: Hybridization expansion results for the current through an interacting dot with U/T = 8, ed + U/2 — 0, and a hard 
bandcutoff uJc/T = 10 (/3r = vT = 10). Left panel: left and right current for V/T — and 5. Il is the current from the 
dot to the left lead and Ir the current from the right lead to the dot. Right panel: average current I = {Il + Ir)/2 and 
dn/dt = Il — Ir for the same parameters. 



state value and the dot occupancy may change. During this transient period, which has been studied in detail in 
Ref. ^H, the current into the dot from the right lead {Ir) need not equal the current out of the dot into the left lead 
{II)- Figure [T51 shows hybridization expansion results for the relaxation dynamics in a dot with voltage bias V/T = 
and 5. As the dot-lead hopping is turned on electrons rush from the leads to the initially empty dot, leading to a fast 
initial rise in the current. For = the current from the dot to the left lead {1^) or the right lead (— /_r) eventually 
vanishes. For > 0, we see that current initially flows into the dot from both sides, but as time is increased II and 
—Ir converge to equal and opposite non-vanishing steady state values. The right hand panel shows the difference 
between the left and right current, which is equal to the derivative of the dot occupation number: Il — Ir = dn/dt, 
with n = -\- ni. This quantity depends relatively weakly on voltage and converges to zero as the steady state is 
reached. 

The average current / — {Il + Ir)/"^ grows with V, as illustrated in the left hand panel of Fig. [111 For the 
parameters in this figure {U/T = 8, €d + U /2 — band cutoff lOc/T — 10, PT ~ vT ~ 10) the small oscillations at 
intermediate times mean that we cannot obtain an accurate estimate of the steady state current. In the right hand 
panel we therefore show the current measured at time tT — 1 (solid lines) and 1.25 (dashed lines) as a function of 
voltage. At large voltage, one observes a slow increase of I{V) in the "Coulomb blockade" regime {V < U) followed 
by a more rapid increase in the current once the voltage bias exceeds the splitting between the Hubbard bands of 
approximately U . We do not find a rapid increase (comparable to the U = curve) in the current near V — Q, 
presumably because the time-scales reached in this simulation are not long enough for a Kondo resonance to form, 
or because the latter is destroyed by even a small applied voltage. However, at voltages V/T > 2, where the Kondo 
resonance is wiped out, we expect our hybridization expansion results to be fairly accurate. 

Figure [T7| compares the hybridization expansion results to the mean field and perturbative calculations. The right 
hand panel shows the comparison to perturbation theory described in Ref. [s^l • We see that for U < 6T the results 
agree quite well. The deviation seen in the U — current is due to a difference in bandwidths {u!c/T = 10 in the Monte 
Carlo simulation, and infinite bandwidth in the analytical calculation). However, at very small V the hybridization 
expansion results indicate a lower current than the perturbation expansion of the self-energy. We believe that this 
difference arises because the hybridization expansion has not been run for long enough times {Tt = 1.25) to capture 
the formation of the Kondo (or fermi liquid) resonance. In the perturbative calculation the crossover from the low V 
un-renormalized behavior to the larger V suppressed I{V) (visible as a flattening of the perturbative I{V) curve at 
V ^ 2T for U — 6T) occurs via a voltage-induced splitting of the Kondo resonance, which was also observed in NCA 
calculations,— However, in these calculations the crossover occurs at a voltage far higher than the Kondo temperature, 
suggesting that the splitting of the Kondo resonance and the associated "hump" in I{V) might be an artifact and 
that further investigation of the crossover would be worthwhile. 

The left panel of Fig. [T7] compares the hybridization expansion and mean field calculations. The data for [/ = 
show that the current measured at = 2 gives a good estimate of the steady state result, especially for larger 
voltage biases (note that the non- interacting model provides a non-trivial test for the strong coupling method). 
The mean field theory clearly underestimates the low V current (and in this regime the QMC data probably are 
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FIG. 16: Hybridization expansion results for the current through an interacting dot with ed + U /2 = G, and a hard bandcutoff 
ujc/T — 10 (PV — i/T — 10). The steady state dot occupancy for this td is 1 (half filling). Left panel: average current 
/ — {II + Ib)/2 for U /T = 8 and indicated values of V . Right hand panel: Current at time tV = 1 (solid lines) and tV = 1.25 
(dashed lines) as a function of voltage, for indicated values of the interaction. The current for the non-interacting dot has been 
measured at = 2. 




FIG. 17: Left panel: comparison to the mean-field result. The circles show the (exact) non-interacting current for lJc/F = 10, 
PV = vF = 10, which is in good agreement with the hybridization expansion result measured at tF = 2, especially for V/T > 3. 
The triangles show the current obtained with the "fixed gap" calculation for U/T = 8, 12 and an effective temperature 
T — 0.05V^ (open symbols) and T = 0.2V (full symbols). Stars show the Monte Carlo results for U/F — 8 and 12 measured at 
tr = 1.25. The right panel compares Monte Carlo results measured at tF = 1 (circles), 1.25 (stars) and 2 (diamonds) to the 
current (for infinite bandwidth) deduced from the 4*'' order perturbation calculation of Ref. [30| . 



themselves an underestimate). However, at larger V the qualitative behavior of the QMC calculations can be more or 
less reproduced by "fixed gap" calculations, if the "effective temperature" (proportional to V) is properly adjusted. 
Figure [T7] also shows that the interacting current approaches the non-interacting value as V becomes very large. The 
comparison provides evidence of the correctness of the simulation results at large biases. It shows in particular that 
we are able to access long enough times to obtain reasonable estimates of the asymptotic behavior and suggests that 
future studies of the pseudothermal broadening effect may be possible. Further investigation of the extent to which 
Coulomb-blockade-like features persist at high bias and strong coupling would also be of interest. 
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VII. CONCLUSIONS 



In this paper we have investigated an approach to nonequilibrium problems based on a stochastic sampling of 
diagrams on the Keldysh contour, with diagrams selected on the basis of their contributions to the expectation value 
being calculated. Both an expansion in interaction strength and an expansion in the dot-lead hybridization were 
considered. The average expansion order scales linearly with the time interval to be studied. In both methods the key 
difficulty is a dynamical sign problem arising from the factors appearing because one expands e^'*^. The average 
sign decays exponentially with perturbation order, and when it becomes smaller than about 0.001 the measurement 
becomes prohibitively difficult. 

For weakly interacting dots, the weak coupling method can be carried to longer times than the hybridization 
expansion. A further advantage of the weak-coupling method is that it starts from an initial density matrix which 
already contains the entanglement between the dot and the leads, so less time is needed to reach steady state. However, 
the growth of the average perturbation order with interaction strength was found to be such that only interactions in 
the weak to intermediate coupling regime U < ttF can be studied. 

The strong coupling method exhibits somewhat worse convergence properties. In interacting dots, only times of 
the order of 1-2 inverse level widths could be reached. A difficulty is that as the method has been formulated here, 
the initial state is a decoupled dot-lead state, which means that the simulation has to build up the necessary dot-lead 
entanglement before steady state can be reached. On the other hand, the method works equally well for all interaction 
strengths and the times accessible appear to be long enough that steady state behavior can be reached, at least at 
large biases where the Kondo effect is not relevant. A general advantage of the diagrammatic Monte Carlo technique 
compared to other methods is that the results are (within the given error bars) exact. There are no discretizations, 
truncations or other approximations. Our results demonstrate that the methods have potential for the simulation of 
more realistic situations, and may also be able to provide basic insights into issues including the crossover from the 
Kondo (unrenormalized differential conductance) to high bias regime. 

We have not attempted to optimize either of the methods. Better choices of cutoff and of initial conditions are 
likely to improve the performance of the algorithms. Better sampling procedures, improved estimators or blocking 
techniques should help reduce the sign problem. Most promising in our opinion are strategies to reduce the average 
perturbation order, for example through the explicit treatment of bath states in the hybridization expansion approach. 
Starting the real-time evolution from a thermalized state by sampling configurations on an "L-shaped" contour with an 
additional branch along the imaginary time direction may lead to a more rapid convergence into the non-equilibrium 
steady state. Efforts in these directions are under way. 
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We rewrite here for convenience the Hamiltonian for a level coupled to two leads, a = L,R (absorbing the Hartree 
shift Un/2 into the definition of the level energy e^) 
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VIII. APPENDIX A: WEAK COUPLING FORMALISM 



A. 



Calculation of d — d Green's function 




(67) 



k,a.a — L,R 



,a.a — L,R 



The physics associated with the couphng to the leads can be reconstructed from 




(68) 



k 



and 




(69) 
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with V the principal value symbol. In the infinite bandwidth, constant density of states limit is constant and 
^a(a;) = 0. 

The coupling to the leads provides a self energy A to the dot Green's function. Because the leads are infinite the 
self energy may be computed in terms of the correlators Gcond of the c electrons with hybridization V ~ 0.^^ In the 
Larkin basis the calculation follows the same lines as the equilibrium^'' one with 

k 

where {5 is the usual positive infinitesimal and the upper (lower) sign pertains to {G^)) 

Gfond,a (fc, ^) = -27r*<5 [Lu^Sk) tanh (^^^f^) ■ (72) 
Inserting Eqs. ([72]) into Eq. jTOl) gives 

AfH = 5„(w)-ir„H, (73) 

^i{Lo) = 5„(c^)+zr„(u;), (74) 
A^(c.) = -2zr„(c.)tanh(^ii^). (75) 

Then using the symbol without the a subscript to denote the sum of left and right channel contributions (so e.g. 
^R,A,K _ ^R:A,K _|_ ^R^A,K -j ^-^jj^ ^ Green's function Gdd is given by 



G'^a G^A -i i ^-^d ^ _ f A^(^) A^(c.) 



(76) 



Use of the standard relationssa. gives 



r,> _ 1 (rK >nR nA\ ^ST -^^^(1 + tanh{{uj - fig,)/ {2Ta,))) 

'-'dd - 2^ dd + ^dd ^dd) - 2^ _ g 5)2 _j_ p2 ' 

^< _ 1 (r^K r^R ,r<A\_sr »ra(l - tanh((Lj - / {2Ta,))) 

^dd — 2 y d.d ~ ^dd + ^dd) — 2^ {uj - Ed - s)'^ + r2 ' ^ ' 



B. calculation of A{t, t') 

We express the quantity A{t,t') = {a\{t')d{t))Q (with retarded/advanced/Keldysh nature here left unspecified) as 

^ = -^Y.^k^Gt■ (79) 

k 

Use of the equation of motion that led to Eq. (11) of Ref. 26] gives (denoting convolution by products) 

\^ A^] [ GiJ\ At) [ Gi,At J ^ ' 

From Eqs. ^ and ^ we find 

_ juj-Ed- S -iV) {Sl -iTL) 

^ " (c. - Ed - sr + n ' (^1) 

aA _ . {i^ - ed - S + iT) (SL+iTL) 

^ - (C. - - 5)2 + n ' (^2) 

.if _ .{u;~ed-S~iT) (-2^^i/^L) - 2z(rL/iL + r^^/i^) (5l + zE^) 

-ri — — 2 



(t^ - Ed - 5)2 + r2 

(-2r^fe^) (c^ - Ed) + 2TLhL {Sr + iPfl) - 2rflfej, (5l + iTl) 

[LO-Ed- 5)2 + r2 



(83) 
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with ha = tanh((w — iia) / i'2Ta)) and T = + F^. Therefore, with /„ = /(cj — fia) denoting the Fermi function for 
lead a and use of the relation A"^ = \{A^ — + A"^) we find 



^ - ^ i.-e,-Sr+T^ ■ 



Setting 5 = gives 



.< _ ^ ^FlFa (/l - /r) ~ TlIl {i^ - £d) f„^. 



C. Mean Field Theory 

In the mean field theory^^ of the nonequilibrium Anderson model one replaces the Hamiltonian by 



(86) 



with 

^o- = £o + f/n_o.. 

The occupancy of the d-orbital of spin a is then 

duj TL{uj)f{uj - /Xi) + TRf{uj - iir) 



(87) 



(^-£, -5(c^))2+r(c^)2 



and one requires self-consistency between Eqs. (1871) and (|88|) . In practice self consistency is achieved by starting from 
an initial guess and iterating until the equations cease to change. 

In our explicit calculations we took a flat band with a hard cutoff defined by 



FL,i?(w) 
Sl,r 



^tan"^ 




+ tan ^ 




) 


S 


S 











O.SF 



In 



(wc — uj)'^ + 6' 
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(89) 
(90) 



with ojc = lOF and 6 = O.IF. 

We choose conventions such that fJ^L + f-^R — and ^ Mfi = V- Then the current is given by 



E 



2rL(tj)r_R(w) (/(w - ^l) - /(t^ - m-r)) 



27r 



(c^-£, -5(co))2+r(c^)2 



(91) 



To represent the broadening effect of a voltage bias, in some calculations we include in the Fermi functions an effective 
temperature equal to a constant times the voltage bias. 
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